{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Yield Management \n",
    "\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "See how mathematical optimization can make your revenues and profits soar in this example, where we’ll show you how an airline can use the AI technology to devise an optimal seat pricing strategy. You’ll learn how to formulate this Yield Management Problem as a three-period stochastic programming problem using the Gurobi Python API and solve it with the Gurobi Optimizer.\n",
    "\n",
    "This model is example 24 from the fifth edition of Model Building in Mathematical Programming by H. Paul Williams on pages 282-284 and 337-340.\n",
    "\n",
    "This modeling example is at the advanced level, where we assume that you know Python and the Gurobi Python API and that you have advanced knowledge of building mathematical optimization models. Typically, the objective function and/or constraints of these examples are complex or require advanced features of the Gurobi Python API.\n",
    "\n",
    "**Download the Repository** <br /> \n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). \n",
    "\n",
    "**Gurobi License** <br /> \n",
    "In order to run this Jupyter Notebook properly, you must have a Gurobi license. If you do not have one, you can request an [evaluation license](https://www.gurobi.com/downloads/request-an-evaluation-license/) as a *commercial user*, or download a [free license](https://www.gurobi.com/academia/academic-program-and-licenses/) as an *academic user*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Description\n",
    "\n",
    "An airline is selling tickets for flights to a particular destination. The flights will depart in three weeks and require up to six planes, each costing £50 000 to hire. Each plane has the following:\n",
    "\n",
    "* 37 First Class seats\n",
    "* 38 Business Class seats\n",
    "* 47 Economy Class seats.\n",
    "\n",
    "The airline needs to decide on an initial price for each of these seats, and will then have the opportunity to update the price after one week and two weeks. Once a customer has purchased a ticket, there is no cancellation option. For administrative simplicity, three price level options are possible in each class (one of which must be chosen). The same option need not be chosen for each class. These are given in the following table for the current period (period 1) and two future periods.\n",
    "\n",
    "![optionsClass](optionsClass.PNG)\n",
    "\n",
    "Demand is uncertain but will be affected by price. Demand forecasts have been made according to a probability distribution that divides the demand levels into three scenarios for each period. The probabilities of the three scenarios\n",
    "in each period are as follows:\n",
    "\n",
    "![scenariosProb](scenariosProb.PNG)\n",
    "\n",
    "The forecasted demand levels are shown in the following table:\n",
    "\n",
    "![forecastDemand](forecastDemand.PNG)\n",
    "\n",
    "The goal is to determine the price levels for the current period, how many seats to sell in each class, the provisional number of planes to book and provisional price levels and seats to sell in future periods in order to maximize\n",
    "expected yield. We should be able to meet commitments under all possible combinations of scenarios. With hindsight (i.e. not known until the beginning of the next period), it turned out that demand in each period (depending on the price level you chose) was as shown in following table.\n",
    "\n",
    "![actualDemand](actualDemand.PNG)\n",
    "\n",
    "We use the actual demand that resulted from the prices we set in period 1 to rerun the model at the beginning of period 2 to set price levels for period 2 and provisional price levels for period 3. We repeat this procedure with a rerun at the beginning of period 3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Model Formulation\n",
    "\n",
    "The Yield Management problem is formulated as a three-period stochastic programming problem. Solving this model for the first time gives recommendations for price levels and sales of week 1 and recommends price levels and sales for subsequent weeks under all possible scenarios. The probabilities of these scenarios will be taken into account in order to maximize expected yield. A week later the model will be rerun, taking into account the committed sales and revenue in the first week, to redetermine recommended prices and sales for the second week (i.e. with ‘recourse’) and the third week under all possible scenarios. The procedure will be repeated again a week later.\n",
    "\n",
    "### Sets and indices\n",
    "\n",
    "$i,j,k \\in \\text{Scenarios}$: Indices and set of scenarios.\n",
    "\n",
    "$h \\in \\text{Options}$: Index and set of price options.\n",
    "\n",
    "$c \\in \\text{Class}$: Index and set of seats categories.\n",
    "\n",
    "### Parameters\n",
    "\n",
    "$\\text{price1}_{c,h} \\in \\mathbb{R}^+$: Price of option $h$ chosen for class $c$ in week 1.\n",
    "\n",
    "$\\text{price2}_{i,c,h} \\in \\mathbb{R}^+$: Price of option $h$ chosen for class $c$ in week 2 as a result of scenario $i$ in week 1.\n",
    "\n",
    "$\\text{price3}_{i,j,c,h} \\in \\mathbb{R}^+$: Price of option $h$ chosen for class $c$ in week 3 as a result of scenario $i$ in week 1, and scenario $j$ in week 2.\n",
    "\n",
    "$\\text{forecast1}_{i,c,h} \\in \\mathbb{R}^+$: Forecast demand in week 1 for class $c$ under price option $h$ and scenario $i$.\n",
    "\n",
    "$\\text{forecast2}_{i,j,c,h} \\in \\mathbb{R}^+$: Forecast demand in week 2 for class $c$ under price option $h$ if scenario $i$ holds in week 1, and scenario $j$ in week 2.\n",
    "\n",
    "$\\text{forecast3}_{i,j,k,c,h} \\in \\mathbb{R}^+$: Forecast demand in week 3 for class $c$ under price option $h$ if  scenario $i$ holds in week 1, scenario $j$ in week 2, and scenario $k$ in week 3.\n",
    "\n",
    "$\\text{prob}_i \\in [0,1]$: Probability of scenario $i$.\n",
    "\n",
    "$\\text{cap}_c \\in \\mathbb{N}$: Capacity per plane for class $c$.\n",
    "\n",
    "$\\text{cost} \\in \\mathbb{R}^+$: Cost to hire a plane.\n",
    "\n",
    "\n",
    "### Decision Variables\n",
    "\n",
    "$p1_{c,h} \\in \\{0, 1\\}$: This binary variable is equal to one if price of option $h$ is chosen for class $c$ in week 1.\n",
    "\n",
    "$p2_{i,c,h} \\in \\{0, 1\\}$: This binary variable is equal to one if price of option $h$ is chosen for class $c$ in week 2 as a result of scenario $i$ in week 1.\n",
    "\n",
    "$p3_{i,j,c,h} \\in \\{0, 1\\}$: This binary variable is equal to one if price of option $h$ is chosen for class $c$ in week 3 as a result of scenario $i$ in week 1 and scenario $j$ in week 2.\n",
    "\n",
    "$s1_{i,c,h} \\in \\mathbb{R}^+$: Number of tickets to be sold in week 1 for class $c$ under price option $h$ and scenario $i$.\n",
    "\n",
    "$s2_{i,j,c,h} \\in \\mathbb{R}^+$: Number of tickets to be sold in week 2 for class $c$ under price option $h$ if scenario $i$ holds in week 1, and scenario $j$ in week 2.\n",
    "\n",
    "$s3_{i,j,k,c,h} \\in \\mathbb{R}^+$: Number of tickets to be sold in week 3 for class $c$ under price option $h$ if  scenario $i$ in week 1, scenario $j$ in week 2, and scenario $k$ in week 3.\n",
    "\n",
    "$n \\in \\mathbb{N}$: Number of planes to fly.\n",
    "\n",
    "### Objective function\n",
    "\n",
    "**Profit**: Maximize expected profit.\n",
    "\n",
    "$$\n",
    "\\text{Maximize} \\quad \\text{profit} = (\n",
    "\\sum_{i \\in \\text{Scenarios}} \\sum_{c \\in \\text{Class}} \n",
    "\\sum_{h \\in \\text{Options}}{\\text{prob}_i * \\text{price1}_{c,h} * p1_{c,h} * s1_{i,c,h}  }  +\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\sum_{i \\in \\text{Scenarios}} \\sum_{j \\in \\text{Scenarios}} \\sum_{c \\in \\text{Class}} \\sum_{h \\in \\text{Options}}\n",
    "{\\text{prob}_i *\\text{prob}_j * \\text{price2}_{c,h} * p2_{i,c,h} * s2_{i,j,c,h} } +\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\sum_{i \\in \\text{Scenarios}} \\sum_{j \\in \\text{Scenarios}} \\sum_{k \\in \\text{Scenarios}} \n",
    "\\sum_{c \\in \\text{Class}} \\sum_{h \\in \\text{Options}}\n",
    "{\\text{prob}_i * \\text{prob}_j * \\text{prob}_k * \\text{price3}_{c,h} * p2_{i,j,c,h} * s3_{i,j,k,c,h}} )\n",
    "- \\text{cost} * n\n",
    "$$\n",
    "\n",
    "### Constraints\n",
    "\n",
    "**Price option week 1**: Only one price option must be chosen in each class in week 1.\n",
    "\n",
    "$$\n",
    "\\sum_{h \\in \\text{Options}} p1_{c,h} = 1 \\quad \\forall c \\in \\text{Class}\n",
    "$$\n",
    "\n",
    "**Sales week 1**: Sales cannot exceed forecasted demand in week 1.\n",
    "\n",
    "$$\n",
    "s1_{i,c,h} \\leq \\text{forecast1}_{i,c,h} * p1_{c,h},\n",
    "\\quad \\forall i \\in \\text{Scenarios}, \\; c \\in \\text{Class}, \\; h \\in \\text{Options}\n",
    "$$\n",
    "\n",
    "**Price option week 2**: Only one price option must be chosen in each class in week 2 for each scenario in week 1.\n",
    "\n",
    "$$\n",
    "\\sum_{h \\in \\text{Options}} p2_{i,c,h} = 1 \\quad \\forall c \\in \\text{Class}, \\; i \\in \\text{Scenarios}\n",
    "$$\n",
    "\n",
    "**Sales week 2**: Sales cannot exceed forecasted demand in week 2.\n",
    "\n",
    "$$\n",
    "s2_{i,j,c,h} \\leq \\text{forecast2}_{j,c,h} * p2_{i,c,h},\n",
    "\\quad \\forall i,j \\in \\text{Scenarios}, \\; c \\in \\text{Class}, \\; h \\in \\text{Options}\n",
    "$$\n",
    "\n",
    "**Price option week 3**: Only one price option must be chosen in each class in week 3 for each scenario in week 1 and week 2.\n",
    "\n",
    "$$\n",
    "\\sum_{h \\in \\text{Options}} p3_{i,j,c,h} = 1 \\quad \\forall c \\in \\text{Class}, \\; i,j \\in \\text{Scenarios}\n",
    "$$\n",
    "\n",
    "**Sales week 3**: Sales cannot exceed forecasted demand in week 3.\n",
    "$$\n",
    "s3_{i,j,k,c,h} \\leq \\text{forecast3}_{k,c,h} * p3_{i,j,c,h},\n",
    "\\quad \\forall i,j,k \\in \\text{Scenarios}, \\; c \\in \\text{Class}, \\; h \\in \\text{Options}\n",
    "$$\n",
    "\n",
    "**Class capacity**: Capacity constraint for each class.\n",
    "\n",
    "$$\n",
    "\\sum_{h \\in \\text{Options}} s1_{i,c,h} +\n",
    "\\sum_{h \\in \\text{Options}} s2_{i,j,c,h} + \n",
    "\\sum_{h \\in \\text{Options}} s3_{i,j,k,c,h} \\leq \\text{cap}_c * n\n",
    "\\quad \\forall i,j,k \\in \\text{Scenarios}, \\; c \\in \\text{Class}\n",
    "$$\n",
    "\n",
    "**Planes** Up to six planes can be hired.\n",
    "\n",
    "$$\n",
    "n \\leq 6\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Python Implementation\n",
    "\n",
    "We import the Gurobi Python Module and other Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy\n",
    "# Note, the restricted license is not sufficient to run this notebook, a full license is needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.7.0 & Gurobi 9.1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input Data\n",
    "We define all the input data for the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lists of classes, price options, and scenarios.\n",
    "\n",
    "classes = ['First', 'Business', 'Economy']\n",
    "\n",
    "options = ['option1', 'option2', 'option3']\n",
    "\n",
    "scenarios = ['sce1', 'sce2', 'sce3']\n",
    "\n",
    "#  Classes, price options, and prices value for week 1, 2, and 3.\n",
    "\n",
    "ch, price1, price2, price3  = gp.multidict({\n",
    "    ('First', 'option1'): [1200, 1400, 1500],\n",
    "    ('Business', 'option1'): [900, 1100, 820],\n",
    "    ('Economy', 'option1'): [500, 700, 480],\n",
    "    ('First', 'option2'): [1000, 1300, 900],\n",
    "    ('Business', 'option2'): [800, 900, 800],\n",
    "    ('Economy', 'option2'): [300, 400, 470],\n",
    "    ('First', 'option3'): [950, 1150, 850],\n",
    "    ('Business', 'option3'): [600, 750, 500],\n",
    "    ('Economy', 'option3'): [200, 350, 450]\n",
    "})\n",
    "\n",
    "# Probablity of each scenario\n",
    "\n",
    "prob ={'sce1': 0.1, 'sce2': 0.7, 'sce3': 0.2}\n",
    "\n",
    "# Forecasted demand for each class, price option, and scenario at week 1\n",
    "\n",
    "ich, fcst1, fcst2, fcst3 = gp.multidict({\n",
    "    ('sce1', 'First', 'option1'): [10, 20, 30],\n",
    "    ('sce1', 'Business', 'option1'): [20, 42, 40],\n",
    "    ('sce1', 'Economy', 'option1'): [45, 50, 50],\n",
    "    ('sce1', 'First', 'option2'): [15, 25, 35],\n",
    "    ('sce1', 'Business', 'option2'): [25, 45, 50],\n",
    "    ('sce1', 'Economy', 'option2'): [55, 52, 60],\n",
    "    ('sce1', 'First', 'option3'): [20, 35, 40],\n",
    "    ('sce1', 'Business', 'option3'): [35, 46, 55],\n",
    "    ('sce1', 'Economy', 'option3'): [60, 60, 80],\n",
    "    ('sce2', 'First', 'option1'): [20, 10, 30],\n",
    "    ('sce2', 'Business', 'option1'): [40, 50, 10],\n",
    "    ('sce2', 'Economy', 'option1'): [50, 60, 50],\n",
    "    ('sce2', 'First', 'option2'): [25, 40, 40],\n",
    "    ('sce2', 'Business', 'option2'): [42, 60, 40],\n",
    "    ('sce2', 'Economy', 'option2'): [52, 65, 60],\n",
    "    ('sce2', 'First', 'option3'): [35, 50, 60],\n",
    "    ('sce2', 'Business', 'option3'): [45, 80, 45],\n",
    "    ('sce2', 'Economy', 'option3'): [63, 90, 70],\n",
    "    ('sce3', 'First', 'option1'): [45, 50, 50],\n",
    "    ('sce3', 'Business', 'option1'): [45, 20, 40],\n",
    "    ('sce3', 'Economy', 'option1'): [55, 10, 60],\n",
    "    ('sce3', 'First', 'option2'): [50, 55, 70],\n",
    "    ('sce3', 'Business', 'option2'): [46, 30, 45],\n",
    "    ('sce3', 'Economy', 'option2'): [56, 40, 65],\n",
    "    ('sce3', 'First', 'option3'): [60, 80, 80],\n",
    "    ('sce3', 'Business', 'option3'): [47, 50, 60],\n",
    "    ('sce3', 'Economy', 'option3'): [64, 60, 70]    \n",
    "})\n",
    "\n",
    "#  Actual demand at weeks 1, 2, and 3.\n",
    "\n",
    "ch, demand1, demand2, demand3  = gp.multidict({\n",
    "    ('First', 'option1'): [25, 22, 45],\n",
    "    ('Business', 'option1'): [50, 45, 20],\n",
    "    ('Economy', 'option1'): [50, 50, 55],\n",
    "    ('First', 'option2'): [30, 45, 60],\n",
    "    ('Business', 'option2'): [40, 55, 40],\n",
    "    ('Economy', 'option2'): [53, 60, 60],\n",
    "    ('First', 'option3'): [40, 50, 75],\n",
    "    ('Business', 'option3'): [45, 75, 50],\n",
    "    ('Economy', 'option3'): [65, 80, 75]\n",
    "})\n",
    "\n",
    "# Class capacity\n",
    "\n",
    "cap ={'First': 37, 'Business': 38, 'Economy': 47}\n",
    "\n",
    "# cost per plane\n",
    "\n",
    "cost = 50000\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preprocessing\n",
    "\n",
    "# week 1 data structures\n",
    "list_ch = []\n",
    "\n",
    "for c in classes:\n",
    "    for h in options:\n",
    "        tp = c,h\n",
    "        list_ch.append(tp)\n",
    "\n",
    "ch = gp.tuplelist(list_ch)\n",
    "\n",
    "list_ich = []\n",
    "\n",
    "for i in scenarios:\n",
    "    for c in classes:\n",
    "        for h in options:\n",
    "            tp = i,c,h\n",
    "            list_ich.append(tp)\n",
    "\n",
    "ich = gp.tuplelist(list_ich)\n",
    "\n",
    "# week 2 data structure\n",
    "\n",
    "list_ijch = []\n",
    "\n",
    "for i in scenarios:\n",
    "    for j in scenarios:\n",
    "        for c in classes:\n",
    "            for h in options:\n",
    "                tp = i,j,c,h\n",
    "                list_ijch.append(tp)\n",
    "\n",
    "ijch = gp.tuplelist(list_ijch)\n",
    "\n",
    "# week 3 data structure\n",
    "\n",
    "list_ijkch = []\n",
    "\n",
    "for i in scenarios:\n",
    "    for j in scenarios:\n",
    "        for k in scenarios:\n",
    "            for c in classes:\n",
    "                for h in options:\n",
    "                    tp = i,j,k,c,h\n",
    "                    list_ijkch.append(tp)\n",
    "\n",
    "ijkch = gp.tuplelist(list_ijkch)\n",
    "\n",
    "# scenarios data structure\n",
    "\n",
    "list_ijk = []\n",
    "\n",
    "for i in scenarios:\n",
    "    for j in scenarios:\n",
    "        for k in scenarios:\n",
    "            tp = i,j,k,\n",
    "            list_ijk.append(tp)\n",
    "\n",
    "ijk = gp.tuplelist(list_ijk)\n",
    "\n",
    "# capacity constraints data structure\n",
    "\n",
    "list_ijkc = []\n",
    "\n",
    "for i in scenarios:\n",
    "    for j in scenarios:\n",
    "        for k in scenarios:\n",
    "            for c in classes:\n",
    "                tp = i,j,k,c\n",
    "                list_ijkc.append(tp)\n",
    "\n",
    "ijkc = gp.tuplelist(list_ijkc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Deployment\n",
    "\n",
    "Solving bilinear problems with Gurobi is as easy as configuring the global parameter `nonConvex`, and setting this parameter to the value of 2.\n",
    "\n",
    "### First Period Model\n",
    "\n",
    "At the beginning of week 1, we want to determine the price options for this week.\n",
    "\n",
    "We create a model and the variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi\\gurobi.lic\n",
      "Changed value of parameter nonConvex to 2\n",
      "   Prev: -1  Min: -1  Max: 2  Default: -1\n"
     ]
    }
   ],
   "source": [
    "model = gp.Model('YieldManagement')\n",
    "\n",
    "# Set global parameters. \n",
    "model.params.nonConvex = 2\n",
    "\n",
    "# Decision variables\n",
    "\n",
    "# price option binary variables at each week\n",
    "p1ch = model.addVars(ch, vtype=GRB.BINARY, name=\"p1ch\")\n",
    "p2ich = model.addVars(ich, vtype=GRB.BINARY, name=\"p2ich\")\n",
    "p3ijch = model.addVars(ijch, vtype=GRB.BINARY, name=\"p3ijch\")\n",
    "\n",
    "# tickets to be sold at each week\n",
    "s1ich = model.addVars(ich, name=\"s1ich\")\n",
    "s2ijch = model.addVars(ijch, name=\"s2ijch\")\n",
    "s3ijkch = model.addVars(ijkch, name=\"s3ijkch\")\n",
    "\n",
    "# number of planes to fly\n",
    "n = model.addVar(ub=6, vtype=GRB.INTEGER, name=\"planes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Week 1 constraints\n",
    "\n",
    "The following constraints ensure that only one price option is chosen in each class in week 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Price option constraints for week 1\n",
    "\n",
    "priceOption1 = model.addConstrs( (gp.quicksum(p1ch[c,h] for h in options ) == 1 for c in classes ), name='priceOption1' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following constraints enforce that sales cannot exceed forecasted demand in week 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sales constraints for week 1\n",
    "\n",
    "sales1 = model.addConstrs( (s1ich[i,c,h] <= fcst1[i,c,h]*p1ch[c,h] for i,c,h in ich ), name='sales1' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Week 2 constraints\n",
    "\n",
    "Only one price option must be chosen in each class in week 2 for each scenario in week 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Price option constraints for week 2\n",
    "\n",
    "priceOption2 = model.addConstrs( (gp.quicksum(p2ich[i,c,h] for h in options ) \n",
    "                                  == 1 for i in scenarios for c in classes ), name='priceOption2' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sales cannot exceed forecasted demand in week 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sales constraints for week 2\n",
    "\n",
    "sales2 = model.addConstrs( (s2ijch[i,j,c,h] <= fcst2[j,c,h]*p2ich[i,c,h] for i,j,c,h in ijch ), name='sales2' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Week 3 constraints\n",
    "\n",
    "Only one price option must be chosen in each class in week 3 for each scenario in week 1 and week 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Price option constraints for week 3\n",
    "\n",
    "priceOption3 = model.addConstrs( (gp.quicksum(p3ijch[i,j,c,h] for h in options ) \n",
    "                                  == 1 for i in scenarios for j in scenarios for c in classes ), name='priceOption3' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sales cannot exceed forecasted demand in week 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sales constraints for week 3\n",
    "\n",
    "sales3 = model.addConstrs( (s3ijkch[i,j,k,c,h] <= fcst3[k,c,h]*p3ijch[i,j,c,h] for i,j,k,c,h in ijkch ), name='sales3' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Capacity constraint for each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Class capacity constraints\n",
    "\n",
    "classCap = model.addConstrs( (gp.quicksum(s1ich[i,c,h] for h in options)  \n",
    "                              + gp.quicksum(s2ijch[i,j,c,h] for h in options) \n",
    "                              + gp.quicksum(s3ijkch[i,j,k,c,h] for h in options)  \n",
    "                              <= cap[c]*n for i,j,k,c in ijkc ) , name='classCap')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective is to maximize expected profit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Objective function\n",
    "obj = gp.quicksum(prob[i]*price1[c,h]*p1ch[c,h]*s1ich[i,c,h] for i,c,h in ich ) \\\n",
    "+ gp.quicksum(prob[i]*prob[j]*price2[c,h]*p2ich[i,c,h]*s2ijch[i,j,c,h] for i,j,c,h in ijch ) \\\n",
    "+ gp.quicksum(prob[i]*prob[j]*prob[k]*price3[c,h]*p3ijch[i,j,c,h]*s3ijkch[i,j,k,c,h] for i,j,k,c,h in ijkch) - cost*n\n",
    "\n",
    "model.setObjective( obj, GRB.MAXIMIZE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 471 rows, 469 columns and 1629 nonzeros\n",
      "Model fingerprint: 0x21b1ab1b\n",
      "Model has 351 quadratic objective terms\n",
      "Variable types: 351 continuous, 118 integer (117 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 9e+01]\n",
      "  Objective range  [5e+04, 5e+04]\n",
      "  QObjective range [9e-01, 2e+03]\n",
      "  Bounds range     [1e+00, 6e+00]\n",
      "  RHS range        [1e+00, 1e+00]\n",
      "Found heuristic solution: objective -300000.0000\n",
      "Presolve time: 0.00s\n",
      "Presolved: 822 rows, 820 columns, 2682 nonzeros\n",
      "Variable types: 702 continuous, 118 integer (117 binary)\n",
      "\n",
      "Root relaxation: objective -1.704970e+05, 487 iterations, 0.01 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 170497.004    0   15 -300000.00 170497.004   157%     -    0s\n",
      "H    0     0                    150842.50000 170497.004  13.0%     -    0s\n",
      "H    0     0                    163722.70000 170497.004  4.14%     -    0s\n",
      "     0     0 168398.617    0   32 163722.700 168398.617  2.86%     -    0s\n",
      "H    0     0                    164667.90000 168398.617  2.27%     -    0s\n",
      "     0     0 166391.907    0   29 164667.900 166391.907  1.05%     -    0s\n",
      "     0     0 166355.400    0   29 164667.900 166355.400  1.02%     -    0s\n",
      "H    0     0                    166150.82456 166355.400  0.12%     -    0s\n",
      "     0     0 166215.053    0   33 166150.825 166215.053  0.04%     -    0s\n",
      "H    0     0                    166195.80000 166215.053  0.01%     -    0s\n",
      "     0     0 166208.253    0   27 166195.800 166208.253  0.01%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 8\n",
      "  MIR: 50\n",
      "  Flow cover: 1\n",
      "\n",
      "Explored 1 nodes (705 simplex iterations) in 0.10 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 6: 166196 166151 164668 ... -300000\n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 1.661958000000e+05, best bound 1.662082533333e+05, gap 0.0075%\n",
      "\n",
      "\n",
      "\n",
      "____________________Week 1 solution___________________________\n",
      "The expected total profit is: £ 166,195.8\n",
      "Number of planes to book: 3.0\n",
      "\n",
      "____________________Week 1 prices_______________________________\n",
      "(First,option1) = £ 1,200\n",
      "(Business,option1) = £ 900\n",
      "(Economy,option1) = £ 500\n",
      "\n",
      "_____________Week 2 provisional prices____________________________\n",
      "(sce1, First, option3) = £ 1,150\n",
      "(sce1, Business, option1) = £ 1,100\n",
      "(sce1, Economy, option1) = £ 700\n",
      "(sce2, First, option3) = £ 1,150\n",
      "(sce2, Business, option1) = £ 1,100\n",
      "(sce2, Economy, option1) = £ 700\n",
      "(sce3, First, option2) = £ 1,300\n",
      "(sce3, Business, option1) = £ 1,100\n",
      "(sce3, Economy, option1) = £ 700\n",
      "\n",
      "_____________Week 3 provisional prices____________________________\n",
      "(sce1, sce1, First, option1) = £ 1,500\n",
      "(sce1, sce1, Business, option2) = £ 800\n",
      "(sce1, sce1, Economy, option1) = £ 480\n",
      "(sce1, sce2, First, option1) = £ 1,500\n",
      "(sce1, sce2, Business, option2) = £ 800\n",
      "(sce1, sce2, Economy, option1) = £ 480\n",
      "(sce1, sce3, First, option1) = £ 1,500\n",
      "(sce1, sce3, Business, option2) = £ 800\n",
      "(sce1, sce3, Economy, option3) = £ 450\n",
      "(sce2, sce1, First, option1) = £ 1,500\n",
      "(sce2, sce1, Business, option2) = £ 800\n",
      "(sce2, sce1, Economy, option1) = £ 480\n",
      "(sce2, sce2, First, option1) = £ 1,500\n",
      "(sce2, sce2, Business, option2) = £ 800\n",
      "(sce2, sce2, Economy, option1) = £ 480\n",
      "(sce2, sce3, First, option1) = £ 1,500\n",
      "(sce2, sce3, Business, option2) = £ 800\n",
      "(sce2, sce3, Economy, option3) = £ 450\n",
      "(sce3, sce1, First, option1) = £ 1,500\n",
      "(sce3, sce1, Business, option2) = £ 800\n",
      "(sce3, sce1, Economy, option1) = £ 480\n",
      "(sce3, sce2, First, option1) = £ 1,500\n",
      "(sce3, sce2, Business, option2) = £ 800\n",
      "(sce3, sce2, Economy, option1) = £ 480\n",
      "(sce3, sce3, First, option1) = £ 1,500\n",
      "(sce3, sce3, Business, option2) = £ 800\n",
      "(sce3, sce3, Economy, option3) = £ 450\n"
     ]
    }
   ],
   "source": [
    "# Verify model formulation\n",
    "\n",
    "model.write('YieldManagement.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model.optimize()\n",
    "\n",
    "#############################################################\n",
    "#            Print results of model for week 1\n",
    "#############################################################\n",
    "\n",
    "print(\"\\n\\n\\n____________________Week 1 solution___________________________\")\n",
    "\n",
    "print(f\"The expected total profit is: £{round(model.objVal,2): ,}\") \n",
    "print(f\"Number of planes to book: {n.x}\")\n",
    "\n",
    "# Week 1 prices for seat class\n",
    "\n",
    "# optimal values of option prices at week 1\n",
    "opt_p1ch = {}\n",
    "\n",
    "print(\"\\n____________________Week 1 prices_______________________________\")\n",
    "for c,h in ch:\n",
    "    opt_p1ch[c,h] = 0\n",
    "    if p1ch[c,h].x > 0.5:\n",
    "        opt_p1ch[c,h] = round(p1ch[c,h].x)\n",
    "        price_ch = opt_p1ch[c,h]*price1[c,h]\n",
    "        print(f\"({c},{h}) = £{price_ch: ,}\")\n",
    "#\n",
    "\n",
    "# Week 2 provisional prices\n",
    "\n",
    "print(\"\\n_____________Week 2 provisional prices____________________________\")\n",
    "for i,c,h in ich:\n",
    "    if p2ich[i,c,h].x > 0.5:\n",
    "        price_ch = round(p2ich[i,c,h].x)*price2[c,h]\n",
    "        print(f\"({i}, {c}, {h}) = £{price_ch: ,}\")\n",
    "\n",
    "# Week 3 provisional prices\n",
    "\n",
    "print(\"\\n_____________Week 3 provisional prices____________________________\")\n",
    "for i,j,c,h in ijch:\n",
    "    if p3ijch[i,j,c,h].x > 0.5:\n",
    "        price_ch = round(p3ijch[i,j,c,h].x)*price3[c,h]\n",
    "        print(f\"({i}, {j}, {c}, {h}) = £{price_ch: ,}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Second Period model\n",
    "\n",
    "We ran the model using the actual demand of week 1 and determine the price options for week 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Changed value of parameter nonConvex to 2\n",
      "   Prev: -1  Min: -1  Max: 2  Default: -1\n",
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 471 rows, 469 columns and 1611 nonzeros\n",
      "Model fingerprint: 0x92c92b4e\n",
      "Model has 351 quadratic objective terms\n",
      "Variable types: 351 continuous, 118 integer (117 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 9e+01]\n",
      "  Objective range  [5e+04, 5e+04]\n",
      "  QObjective range [9e-01, 2e+03]\n",
      "  Bounds range     [1e+00, 6e+00]\n",
      "  RHS range        [1e+00, 1e+00]\n",
      "Found heuristic solution: objective -300000.0000\n",
      "Presolve removed 30 rows and 27 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 765 rows, 766 columns, 2376 nonzeros\n",
      "Variable types: 657 continuous, 109 integer (108 binary)\n",
      "\n",
      "Root relaxation: objective -1.743903e+05, 456 iterations, 0.01 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 174390.259    0   22 -300000.00 174390.259   158%     -    0s\n",
      "H    0     0                    155574.00000 174390.259  12.1%     -    0s\n",
      "H    0     0                    169901.50000 174390.259  2.64%     -    0s\n",
      "     0     0 171102.200    0   30 169901.500 171102.200  0.71%     -    0s\n",
      "     0     0 170934.033    0   30 169901.500 170934.033  0.61%     -    0s\n",
      "H    0     0                    170801.00000 170934.033  0.08%     -    0s\n",
      "     0     0 170848.733    0   33 170801.000 170848.733  0.03%     -    0s\n",
      "     0     0 170813.360    0   33 170801.000 170813.360  0.01%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 1\n",
      "  MIR: 43\n",
      "\n",
      "Explored 1 nodes (632 simplex iterations) in 0.08 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 4: 170801 169901 155574 -300000 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 1.708010000000e+05, best bound 1.708133600000e+05, gap 0.0072%\n",
      "\n",
      "\n",
      "\n",
      "____________________Week 2 solution___________________________\n",
      "The expected total profit at the beginning of week 2 is: £ 170,801.0\n",
      "Number of planes to book: 3.0\n",
      "\n",
      "_____________Week 2 prices____________________________\n",
      "(First,option3) = £ 1,150\n",
      "(Business,option1) = £ 1,100\n",
      "(Economy,option1) = £ 700\n",
      "\n",
      "_____________Week 3 provisional prices____________________________\n",
      "(sce1, sce1, First, option1) = £ 1500\n",
      "(sce1, sce1, Business, option2) = £ 800\n",
      "(sce1, sce1, Economy, option1) = £ 480\n",
      "(sce1, sce2, First, option1) = £ 1500\n",
      "(sce1, sce2, Business, option2) = £ 800\n",
      "(sce1, sce2, Economy, option1) = £ 480\n",
      "(sce1, sce3, First, option1) = £ 1500\n",
      "(sce1, sce3, Business, option2) = £ 800\n",
      "(sce1, sce3, Economy, option3) = £ 450\n",
      "(sce2, sce1, First, option1) = £ 1500\n",
      "(sce2, sce1, Business, option2) = £ 800\n",
      "(sce2, sce1, Economy, option1) = £ 480\n",
      "(sce2, sce2, First, option1) = £ 1500\n",
      "(sce2, sce2, Business, option2) = £ 800\n",
      "(sce2, sce2, Economy, option1) = £ 480\n",
      "(sce2, sce3, First, option1) = £ 1500\n",
      "(sce2, sce3, Business, option2) = £ 800\n",
      "(sce2, sce3, Economy, option3) = £ 450\n",
      "(sce3, sce1, First, option1) = £ 1500\n",
      "(sce3, sce1, Business, option2) = £ 800\n",
      "(sce3, sce1, Economy, option1) = £ 480\n",
      "(sce3, sce2, First, option1) = £ 1500\n",
      "(sce3, sce2, Business, option2) = £ 800\n",
      "(sce3, sce2, Economy, option1) = £ 480\n",
      "(sce3, sce3, First, option1) = £ 1500\n",
      "(sce3, sce3, Business, option2) = £ 800\n",
      "(sce3, sce3, Economy, option3) = £ 450\n"
     ]
    }
   ],
   "source": [
    "model2 = gp.Model('YieldManagement2')\n",
    "\n",
    "# Set global parameters \n",
    "model2.params.nonConvex = 2\n",
    "\n",
    "# Decision variables\n",
    "\n",
    "# price option binary variables at each week\n",
    "p1ch = model2.addVars(ch, vtype=GRB.BINARY, name=\"p1ch\")\n",
    "\n",
    "# Fix price options of week 1\n",
    "for c,h in ch:\n",
    "    p1ch[c,h].lb = opt_p1ch[c,h]\n",
    "\n",
    "p2ich = model2.addVars(ich, vtype=GRB.BINARY, name=\"p2ich\")\n",
    "p3ijch = model2.addVars(ijch, vtype=GRB.BINARY, name=\"p3ijch\")\n",
    "\n",
    "# tickets to be sold at each week\n",
    "s1ich = model2.addVars(ich, name=\"s1ich\")\n",
    "\n",
    "# use hindsight demand of week 1\n",
    "for i,c,h in ich:\n",
    "    fcst1[i,c,h] = 0\n",
    "    fcst1[i,c,h] = demand1[c,h]*opt_p1ch[c,h] \n",
    "\n",
    "s2ijch = model2.addVars(ijch, name=\"s2ijch\")\n",
    "s3ijkch = model2.addVars(ijkch, name=\"s3ijkch\")\n",
    "\n",
    "# number of planes to fly\n",
    "n = model2.addVar(ub=6, vtype=GRB.INTEGER, name=\"planes\")\n",
    "\n",
    "# Price option constraints for week 1\n",
    "\n",
    "priceOption1 = model2.addConstrs( (gp.quicksum(p1ch[c,h] for h in options ) == 1 for c in classes ), name='priceOption1' )\n",
    "\n",
    "# sales constraints for week 1\n",
    "\n",
    "sales1 = model2.addConstrs( (s1ich[i,c,h] <= fcst1[i,c,h]*p1ch[c,h] for i,c,h in ich ), name='sales1' )\n",
    "\n",
    "# Price option constraints for week 2\n",
    "\n",
    "priceOption2 = model2.addConstrs( (gp.quicksum(p2ich[i,c,h] for h in options ) \n",
    "                                  == 1 for i in scenarios for c in classes ), name='priceOption2' )\n",
    "\n",
    "# sales constraints for week 2\n",
    "\n",
    "sales2 = model2.addConstrs( (s2ijch[i,j,c,h] <= fcst2[j,c,h]*p2ich[i,c,h] for i,j,c,h in ijch ), name='sales2' )\n",
    "\n",
    "# Price option constraints for week 3\n",
    "\n",
    "priceOption3 = model2.addConstrs( (gp.quicksum(p3ijch[i,j,c,h] for h in options ) \n",
    "                                  == 1 for i in scenarios for j in scenarios for c in classes ), name='priceOption3' )\n",
    "\n",
    "# sales constraints for week 3\n",
    "\n",
    "sales3 = model2.addConstrs( (s3ijkch[i,j,k,c,h] <= fcst3[k,c,h]*p3ijch[i,j,c,h] for i,j,k,c,h in ijkch ), name='sales3' )\n",
    "\n",
    "# Class capacity constraints.\n",
    "\n",
    "classCap = model2.addConstrs( (gp.quicksum(s1ich[i,c,h] for h in options)  \n",
    "                              + gp.quicksum(s2ijch[i,j,c,h] for h in options) \n",
    "                              + gp.quicksum(s3ijkch[i,j,k,c,h] for h in options)  \n",
    "                              <= cap[c]*n for i,j,k,c in ijkc ) , name='classCap')\n",
    "\n",
    "# Objective function\n",
    "obj = gp.quicksum(prob[i]*price1[c,h]*p1ch[c,h]*s1ich[i,c,h] for i,c,h in ich ) \\\n",
    "+ gp.quicksum(prob[i]*prob[j]*price2[c,h]*p2ich[i,c,h]*s2ijch[i,j,c,h] for i,j,c,h in ijch ) \\\n",
    "+ gp.quicksum(prob[i]*prob[j]*prob[k]*price3[c,h]*p3ijch[i,j,c,h]*s3ijkch[i,j,k,c,h] for i,j,k,c,h in ijkch) - cost*n\n",
    "\n",
    "model2.setObjective( obj, GRB.MAXIMIZE)\n",
    "\n",
    "# Verify model formulation\n",
    "\n",
    "model2.write('YieldManagement2.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model2.optimize()\n",
    "\n",
    "#############################################################\n",
    "#            Print results of model for week 2\n",
    "#############################################################\n",
    "\n",
    "print(\"\\n\\n\\n____________________Week 2 solution___________________________\")\n",
    "\n",
    "print(f\"The expected total profit at the beginning of week 2 is: £ {round(model2.objVal,2):,}\") \n",
    "print(f\"Number of planes to book: {n.x}\")\n",
    "\n",
    "# Week 2 prices\n",
    "\n",
    "# optimal values of option prices at week 1\n",
    "opt_p2ich = {}\n",
    "\n",
    "print(\"\\n_____________Week 2 prices____________________________\")\n",
    "for i,c,h in ich:\n",
    "    opt_p2ich[i,c,h] = 0\n",
    "    if p2ich[i,c,h].x > 0.5:\n",
    "        opt_p2ich[i,c,h] = round(p2ich[i,c,h].x)\n",
    "        price_ch = opt_p2ich[i,c,h]*price2[c,h]\n",
    "        #print(f\"({i},{c},{h}) = £ {price_ch}\")\n",
    "        if i == 'sce1':\n",
    "            print(f\"({c},{h}) = £{price_ch: ,}\")\n",
    "#\n",
    "\n",
    "# Week 3 provisional prices\n",
    "\n",
    "print(\"\\n_____________Week 3 provisional prices____________________________\")\n",
    "for i,j,c,h in ijch:\n",
    "    if p3ijch[i,j,c,h].x > 0.5:\n",
    "        price_ch = round(p3ijch[i,j,c,h].x)*price3[c,h]\n",
    "        print(f\"({i}, {j}, {c}, {h}) = £ {price_ch}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Third Period model\n",
    "We ran the model using the actual demand of weeks 1 and 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Changed value of parameter nonConvex to 2\n",
      "   Prev: -1  Min: -1  Max: 2  Default: -1\n",
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 471 rows, 469 columns and 1557 nonzeros\n",
      "Model fingerprint: 0x9217fd30\n",
      "Model has 351 quadratic objective terms\n",
      "Variable types: 351 continuous, 118 integer (117 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 8e+01]\n",
      "  Objective range  [5e+04, 5e+04]\n",
      "  QObjective range [9e-01, 2e+03]\n",
      "  Bounds range     [1e+00, 6e+00]\n",
      "  RHS range        [1e+00, 1e+00]\n",
      "Found heuristic solution: objective -300000.0000\n",
      "Presolve removed 120 rows and 108 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 594 rows, 604 columns, 1782 nonzeros\n",
      "Variable types: 522 continuous, 82 integer (81 binary)\n",
      "\n",
      "Root relaxation: objective -1.786023e+05, 389 iterations, 0.01 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 178602.270    0   19 -300000.00 178602.270   160%     -    0s\n",
      "H    0     0                    115595.47092 178602.270  54.5%     -    0s\n",
      "H    0     0                    155270.00000 178602.270  15.0%     -    0s\n",
      "H    0     0                    173270.00000 178602.270  3.08%     -    0s\n",
      "     0     0 173910.600    0   37 173270.000 173910.600  0.37%     -    0s\n",
      "H    0     0                    173680.00000 173910.600  0.13%     -    0s\n",
      "     0     0 173695.160    0   30 173680.000 173695.160  0.01%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 1\n",
      "  MIR: 27\n",
      "\n",
      "Explored 1 nodes (627 simplex iterations) in 0.07 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 5: 173680 173270 155270 ... -300000\n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 1.736800000000e+05, best bound 1.736951600000e+05, gap 0.0087%\n",
      "\n",
      "\n",
      "\n",
      "____________________Week 3 solution___________________________\n",
      "The expected total profit is: £ 173,680.0\n",
      "Number of planes to book: 3.0\n",
      "\n",
      "_____________Week 3 prices____________________________\n",
      "(First, option1) = £ 1,500\n",
      "(Business, option2) = £ 800\n",
      "(Economy, option1) = £ 480\n"
     ]
    }
   ],
   "source": [
    "model3 = gp.Model('YieldManagement3')\n",
    "\n",
    "# Set global parameters. \n",
    "model3.params.nonConvex = 2\n",
    "\n",
    "# Decision variables\n",
    "\n",
    "# price option binary variables at each week\n",
    "p1ch = model3.addVars(ch, vtype=GRB.BINARY, name=\"p1ch\")\n",
    "\n",
    "# Fix price options of week 1\n",
    "\n",
    "for c,h in ch:\n",
    "    p1ch[c,h].lb = opt_p1ch[c,h]\n",
    "\n",
    "p2ich = model3.addVars(ich, vtype=GRB.BINARY, name=\"p2ich\")\n",
    "\n",
    "# Fix price options of week 2\n",
    "\n",
    "for i,c,h in ich:\n",
    "    p2ich[i,c,h].lb = opt_p2ich[i,c,h]\n",
    "\n",
    "p3ijch = model3.addVars(ijch, vtype=GRB.BINARY, name=\"p3ijch\")\n",
    "\n",
    "# tickets to be sold at each week\n",
    "s1ich = model3.addVars(ich, name=\"s1ich\")\n",
    "\n",
    "# use hindsight demand of week 1\n",
    "for i,c,h in ich:\n",
    "    fcst1[i,c,h] = 0\n",
    "    fcst1[i,c,h] = demand1[c,h]*opt_p1ch[c,h]\n",
    "\n",
    "s2ijch = model3.addVars(ijch, name=\"s2ijch\")\n",
    "\n",
    "# use hindsight demand of week 2\n",
    "for j,c,h in ich:\n",
    "    fcst2[j,c,h] = 0\n",
    "    fcst2[j,c,h] = demand2[c,h]*opt_p2ich[j,c,h]\n",
    "\n",
    "s3ijkch = model3.addVars(ijkch, name=\"s3ijkch\")\n",
    "\n",
    "# number of planes to fly\n",
    "n = model3.addVar(ub=6, vtype=GRB.INTEGER, name=\"planes\")\n",
    "\n",
    "# Price option constraints for week 1\n",
    "\n",
    "priceOption1 = model3.addConstrs( (gp.quicksum(p1ch[c,h] for h in options ) == 1 for c in classes ), name='priceOption1' )\n",
    "\n",
    "# sales constraints for week 1\n",
    "\n",
    "sales1 = model3.addConstrs( (s1ich[i,c,h] <= fcst1[i,c,h]*p1ch[c,h] for i,c,h in ich ), name='sales1' )\n",
    "\n",
    "# Price option constraints for week 2\n",
    "\n",
    "priceOption2 = model3.addConstrs( (gp.quicksum(p2ich[i,c,h] for h in options ) \n",
    "                                  == 1 for i in scenarios for c in classes ), name='priceOption2' )\n",
    "\n",
    "# sales constraints for week 2\n",
    "\n",
    "sales2 = model3.addConstrs( (s2ijch[i,j,c,h] <= fcst2[j,c,h]*p2ich[i,c,h] for i,j,c,h in ijch ), name='sales2' )\n",
    "\n",
    "# Price option constraints for week 3\n",
    "\n",
    "priceOption3 = model3.addConstrs( (gp.quicksum(p3ijch[i,j,c,h] for h in options ) \n",
    "                                  == 1 for i in scenarios for j in scenarios for c in classes ), name='priceOption3' )\n",
    "\n",
    "# sales constraints for week 3\n",
    "\n",
    "sales3 = model3.addConstrs( (s3ijkch[i,j,k,c,h] <= fcst3[k,c,h]*p3ijch[i,j,c,h] for i,j,k,c,h in ijkch ), name='sales3' )\n",
    "\n",
    "# Class capacity constraints.\n",
    "\n",
    "classCap = model3.addConstrs( (gp.quicksum(s1ich[i,c,h] for h in options)  \n",
    "                              + gp.quicksum(s2ijch[i,j,c,h] for h in options) \n",
    "                              + gp.quicksum(s3ijkch[i,j,k,c,h] for h in options)  \n",
    "                              <= cap[c]*n for i,j,k,c in ijkc ) , name='classCap')\n",
    "\n",
    "# Objective function\n",
    "obj = gp.quicksum(prob[i]*price1[c,h]*p1ch[c,h]*s1ich[i,c,h] for i,c,h in ich ) \\\n",
    "+ gp.quicksum(prob[i]*prob[j]*price2[c,h]*p2ich[i,c,h]*s2ijch[i,j,c,h] for i,j,c,h in ijch ) \\\n",
    "+ gp.quicksum(prob[i]*prob[j]*prob[k]*price3[c,h]*p3ijch[i,j,c,h]*s3ijkch[i,j,k,c,h] for i,j,k,c,h in ijkch) - cost*n\n",
    "\n",
    "model3.setObjective( obj, GRB.MAXIMIZE)\n",
    "\n",
    "# Verify model formulation\n",
    "\n",
    "model3.write('YieldManagement3.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model3.optimize()\n",
    "\n",
    "#############################################################\n",
    "#            Print results of model for week 3\n",
    "#############################################################\n",
    "\n",
    "\n",
    "print(\"\\n\\n\\n____________________Week 3 solution___________________________\")\n",
    "\n",
    "print(f\"The expected total profit is: £ {round(model3.objVal,2):,}\") \n",
    "print(f\"Number of planes to book: {n.x}\")\n",
    "\n",
    "# Week 3  prices\n",
    "\n",
    "# optimal values of option prices at week 3\n",
    "opt_p3ijch = {}\n",
    "\n",
    "\n",
    "print(\"\\n_____________Week 3 prices____________________________\")\n",
    "for i,j,c,h in ijch:\n",
    "    opt_p3ijch[i,j,c,h] = 0\n",
    "    if p3ijch[i,j,c,h].x > 0.5:\n",
    "        opt_p3ijch[i,j,c,h] = round(p3ijch[i,j,c,h].x)\n",
    "        price_ch = opt_p3ijch[i,j,c,h]*price3[c,h]\n",
    "        if i == 'sce1' and j == 'sce1':\n",
    "            print(f\"({c}, {h}) = £{price_ch: ,}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solution at take-off\n",
    "We ran the model using the actual demand of weeks 1, 2, and 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Changed value of parameter nonConvex to 2\n",
      "   Prev: -1  Min: -1  Max: 2  Default: -1\n",
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 471 rows, 469 columns and 1395 nonzeros\n",
      "Model fingerprint: 0x0aff6c8c\n",
      "Model has 351 quadratic objective terms\n",
      "Variable types: 351 continuous, 118 integer (117 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 6e+01]\n",
      "  Objective range  [5e+04, 5e+04]\n",
      "  QObjective range [9e-01, 2e+03]\n",
      "  Bounds range     [1e+00, 6e+00]\n",
      "  RHS range        [1e+00, 1e+00]\n",
      "Found heuristic solution: objective -0.0000000\n",
      "Presolve removed 470 rows and 458 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 1 rows, 11 columns, 11 nonzeros\n",
      "Variable types: 10 continuous, 1 integer (0 binary)\n",
      "\n",
      "Root relaxation: objective 1.952617e+05, 1 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 195261.702    0    1   -0.00000 195261.702      -     -    0s\n",
      "H    0     0                    180419.39045 195261.702  8.23%     -    0s\n",
      "H    0     0                    181705.53191 195261.702  7.46%     -    0s\n",
      "H    0     0                    184030.00000 195261.702  6.10%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 1\n",
      "\n",
      "Explored 1 nodes (1 simplex iterations) in 0.02 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 2: 184030 -0 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 1.840300000000e+05, best bound 1.840300000000e+05, gap 0.0000%\n",
      "\n",
      "\n",
      "\n",
      "____________________Take off solution___________________________\n",
      "The actual total profit is: £184,030.0\n",
      "Number of planes used: 3.0\n",
      "\n",
      "___________Week 1 seats sold and revenue__________________________\n",
      "First class: 25 seats sold at £1,200: revenue £30,000  \n",
      "Business class: 50 seats sold at £900: revenue £45,000  \n",
      "Economy class: 50 seats sold at £500: revenue £25,000  \n",
      "___________Period 2 seats sold and revenue__________________________\n",
      "First class: 41 seats sold at £1,150: revenue £47,150  \n",
      "Business class: 45 seats sold at £1,100: revenue £49,500  \n",
      "Economy class: 50 seats sold at £700: revenue £35,000  \n",
      "___________Period 3 seats sold and revenue__________________________\n",
      "First class: 45 seats sold at £1,500: revenue £67,500  \n",
      "Business class: 19 seats sold at £800: revenue £15,200  \n",
      "Economy class: 41 seats sold at £480: revenue £19,680  \n"
     ]
    }
   ],
   "source": [
    "model4 = gp.Model('YieldManagement4')\n",
    "\n",
    "# Set global parameters \n",
    "model4.params.nonConvex = 2\n",
    "\n",
    "# Decision variables\n",
    "\n",
    "# price option binary variables at each week\n",
    "p1ch = model4.addVars(ch, vtype=GRB.BINARY, name=\"p1ch\")\n",
    "\n",
    "# Fix price options of week 1\n",
    "\n",
    "for c,h in ch:\n",
    "    p1ch[c,h].lb = opt_p1ch[c,h]\n",
    "\n",
    "p2ich = model4.addVars(ich, vtype=GRB.BINARY, name=\"p2ich\")\n",
    "\n",
    "# Fix price options of week 2\n",
    "\n",
    "for i,c,h in ich:\n",
    "    p2ich[i,c,h].lb = opt_p2ich[i,c,h]\n",
    "\n",
    "p3ijch = model4.addVars(ijch, vtype=GRB.BINARY, name=\"p3ijch\")\n",
    "\n",
    "# Capture one scenario where opt_p3ijch[i,j,c,h] = 1\n",
    "opt_p3kch = {}\n",
    "\n",
    "for i,j,c,h in ijch:\n",
    "    p3ijch[i,j,c,h].lb = opt_p3ijch[i,j,c,h]\n",
    "    opt_p3kch[j,c,h] = 0\n",
    "    if opt_p3ijch[i,j,c,h] == 1:\n",
    "        opt_p3kch[j,c,h] = opt_p3ijch[i,j,c,h] \n",
    "\n",
    "# tickets to be sold at each week\n",
    "s1ich = model4.addVars(ich, name=\"s1ich\")\n",
    "\n",
    "# use hindsight demand of week 1\n",
    "for i,c,h in ich:\n",
    "    fcst1[i,c,h] = 0\n",
    "    fcst1[i,c,h] = demand1[c,h]*opt_p1ch[c,h] \n",
    "\n",
    "s2ijch = model4.addVars(ijch, name=\"s2ijch\")\n",
    "\n",
    "# use hindsight demand of week 2\n",
    "for j,c,h in ich:\n",
    "    fcst2[j,c,h] = 0\n",
    "    fcst2[j,c,h] = demand2[c,h]*opt_p2ich[j,c,h]\n",
    "\n",
    "\n",
    "s3ijkch = model4.addVars(ijkch, name=\"s3ijkch\")\n",
    "\n",
    "# use hindsight demand of week 3\n",
    "for k,c,h in ich:\n",
    "    fcst3[k,c,h] = 0\n",
    "    fcst3[k,c,h] = demand3[c,h]*opt_p3kch[k,c,h]\n",
    "\n",
    "    \n",
    "# number of planes to fly\n",
    "n = model4.addVar(ub=6, vtype=GRB.INTEGER, name=\"planes\")\n",
    "\n",
    "# Price option constraints for week 1\n",
    "\n",
    "priceOption1 = model4.addConstrs( (gp.quicksum(p1ch[c,h] for h in options ) == 1 for c in classes ), name='priceOption1' )\n",
    "\n",
    "# sales constraints for week 1\n",
    "\n",
    "sales1 = model4.addConstrs( (s1ich[i,c,h] <= fcst1[i,c,h]*p1ch[c,h] for i,c,h in ich ), name='sales1' )\n",
    "\n",
    "# Price option constraints for week 2\n",
    "\n",
    "priceOption2 = model4.addConstrs( (gp.quicksum(p2ich[i,c,h] for h in options ) \n",
    "                                  == 1 for i in scenarios for c in classes ), name='priceOption2' )\n",
    "\n",
    "# sales constraints for week 2\n",
    "\n",
    "sales2 = model4.addConstrs( (s2ijch[i,j,c,h] <= fcst2[j,c,h]*p2ich[i,c,h] for i,j,c,h in ijch ), name='sales2' )\n",
    "\n",
    "# Price option constraints for week 3\n",
    "\n",
    "priceOption3 = model4.addConstrs( (gp.quicksum(p3ijch[i,j,c,h] for h in options ) \n",
    "                                  == 1 for i in scenarios for j in scenarios for c in classes ), name='priceOption3' )\n",
    "\n",
    "# sales constraints for week 3\n",
    "\n",
    "sales3 = model4.addConstrs( (s3ijkch[i,j,k,c,h] <= fcst3[k,c,h]*p3ijch[i,j,c,h] for i,j,k,c,h in ijkch ), name='sales3' )\n",
    "\n",
    "# Class capacity constraints.\n",
    "\n",
    "classCap = model4.addConstrs( (gp.quicksum(s1ich[i,c,h] for h in options)  \n",
    "                              + gp.quicksum(s2ijch[i,j,c,h] for h in options) \n",
    "                              + gp.quicksum(s3ijkch[i,j,k,c,h] for h in options)  \n",
    "                              <= cap[c]*n for i,j,k,c in ijkc ) , name='classCap')\n",
    "\n",
    "# Objective function\n",
    "obj = gp.quicksum(prob[i]*price1[c,h]*p1ch[c,h]*s1ich[i,c,h] for i,c,h in ich ) \\\n",
    "+ gp.quicksum(prob[i]*prob[j]*price2[c,h]*p2ich[i,c,h]*s2ijch[i,j,c,h] for i,j,c,h in ijch ) \\\n",
    "+ gp.quicksum(prob[i]*prob[j]*prob[k]*price3[c,h]*p3ijch[i,j,c,h]*s3ijkch[i,j,k,c,h] for i,j,k,c,h in ijkch) - cost*n\n",
    "\n",
    "model4.setObjective( obj, GRB.MAXIMIZE)\n",
    "\n",
    "# Verify model formulation\n",
    "\n",
    "model4.write('YieldManagement4.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model4.optimize()\n",
    "\n",
    "#############################################################\n",
    "#            Print results of model for week 4\n",
    "#############################################################\n",
    "\n",
    "\n",
    "print(\"\\n\\n\\n____________________Take off solution___________________________\")\n",
    "\n",
    "print(f\"The actual total profit is: £{round(model4.objVal,2):,}\") \n",
    "print(f\"Number of planes used: {n.x}\")\n",
    "\n",
    "\n",
    "# Sales week 1\n",
    "\n",
    "print(\"\\n___________Week 1 seats sold and revenue__________________________\")\n",
    "for i,c,h in ich:\n",
    "    if i == 'sce1':\n",
    "        if s1ich[i,c,h].x > 1e-6:\n",
    "            tickets = round(s1ich[i,c,h].x)\n",
    "            price = price1[c,h]*round(p1ch[c,h].x)\n",
    "            revenue = price*tickets\n",
    "            print(f\"{c} class: {tickets} seats sold at £{price:,}: revenue £{revenue:,}  \")\n",
    "            \n",
    "# Sales week 2\n",
    "\n",
    "print(\"___________Period 2 seats sold and revenue__________________________\")\n",
    "for i,j,c,h in ijch:\n",
    "    if i == 'sce1' and j == 'sce1':\n",
    "        if s2ijch[i,j,c,h].x > 1e-6:\n",
    "            tickets = round(s2ijch[i,j,c,h].x)\n",
    "            price = price2[c,h]*round(p2ich[i,c,h].x)\n",
    "            revenue = price*tickets\n",
    "            print(f\"{c} class: {tickets} seats sold at £{price:,}: revenue £{revenue:,}  \")\n",
    "            \n",
    "# Sales week 3\n",
    "\n",
    "print(\"___________Period 3 seats sold and revenue__________________________\")\n",
    "for i,j,k,c,h in ijkch:\n",
    "    if i == 'sce1' and j == 'sce1' and k == 'sce1':\n",
    "        if s3ijkch[i,j,k,c,h].x > 1e-6:\n",
    "            tickets = round(s3ijkch[i,j,k,c,h].x )\n",
    "            price = price3[c,h]*round(p3ijch[i,j,c,h].x)\n",
    "            revenue = price*tickets\n",
    "            print(f\"{c} class: {tickets} seats sold at £{price:,}: revenue £{revenue:,}  \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
